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Abstract 

A chemotaxis-diffusion-convection coupling system for describing a form of buoyant convection in 
which the fluid develops convection cells and plume patterns will be investigated numerically in this 
study. Based on the two-dimensional convective chemotaxis-fluid model proposed in the literature, we 
developed an upwind finite element method to investigate the pattern formation and the hydrodynamical 
stability of the system. The numerical simulations illustrate different predicted physical regimes in the 
system. In the convective regime, the predicted plumes resemble Benard instabilities. Our numerical 
results show how structured layers of bacteria are formed before bacterium rich plumes fall in the fluid. 
The plumes have a well defined spectrum of wavelengths and have an exponential growth rate, yet their 
position can only be predicted in very simple examples. In the chemotactic and diffusive regimes, the 
effects of chemotaxis are investigated. Our results indicate that the chemotaxis can stabilize the overall 
system. A time scale analysis has been performed to demonstrate that the critical taxis Rayleigh number 
for which instabilities set in depends on the chemotaxis head and sensitivity. In addition, the comparison 
of the differential systems of chemotaxis-diffusion-convection, double diffusive convection, and Rayleigh- 
Benard convection establishes a set of evidences that even if the physical mechanisms are different at the 
same time the dimensionless systems are strongly related to each other. 
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Table 1: Nomenclature description 


c 

concentration of oxygen 

Greek symbols 


^air 

concentration of oxygen in air 

Ps 

solute expansion coefficient 

D 

diffusivity 

Pt 

thermal expansion coefficient 

G 

growth rate of the plume ampli¬ 

K 

bacterial oxygen consumption 


tude 


rate 

h 

container height 

t 

dimensionless container width 

H 

chemotaxis head 

R 

dynamic viscosity 

j 

vertical unit vector upwards 

V 

kinematic viscosity 

Le 

Lewis number 

P 

fluid density 

n 

areal number density of bacteria 

Pb 

volumetric mass density of a bac¬ 
terium 

n 

unit outward normal vector 



n 0 

initial number density of bacteria 

Superscripts 


no 

initial average number density of 
bacteria 


dimensionless quantities 

V 

pressure 



p r 

Prandtl number 

Subscripts 


Ra 

Rayleigh number 

‘b 

bacterium 

S 

dimensionless chemotaxis sensi¬ 
tivity 

* c 

critical 

S dim 

dimensional chemotaxis sensitiv¬ 
ity 

'conv 

convection 

S 

solute concentration 

Miff 

diffusion 

t 

time 

‘ra 

mass 

T 

time scale of bacterium transport 

‘O 

oxygen 

T 

temperature 

's 

solute 

u 

velocity vector 

*T 

taxis 

u 

velocity in x-direction 



V 

velocity in y-direction 



v b 

volume of a bacterium 



x = (x,y) 

coordinate axes 
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1 Introduction 


Taxis refers to the collective motion of cells or an organism in response to an attractant gradient. The 
nature of the attractant stimulus can be of chemical (chemotaxis), physical (baro-, electro-, magneto-, 
phono-, photo-, and thermotaxis), or mechanical (hapto- and rheotaxis) origins. Chemotaxis refers to cell 
movement primed by a external chemical signal that can either be emitted by the same population of cells 
or created by an external source [I]. In particular, aerotaxis is related to the movement toward a gradient 
of increasing oxygen concentration. 

The phenomenon that couples chemotaxis, diffusion, and convection has been illustrated by experiments 
on suspensions of bacteria in a container filled with water m- Oxygen diffuses in the container from the 
surface. As bacteria consume oxygen, oxygen concentration falls everywhere except at the surface, hence 
creating a vertical concentration difference. Bacteria move up to higher concentration of oxygen and quickly 
get densely packed below the surface in a relatively thin layer. Subsequently, Rayleigh-Benard-like instability 
appears near the surface. These dynamical instabilities exhibit complex convection patterns with plumes of 
bacteria falling in the fluid. 

Chemotaxis-diffusion-convection is a particular case of the so-called bioconvection. Bioconvection is a 
more general term for suspensions of swimming microorganisms which are denser than the solvent fluid. 
Bioconvection and the different mechanisms of upswimming have been reviewed [1]. 

Mathematical modeling of chemotaxis has been introduced by authors in [S| and 0. Mathematical 
analysis mainly focused on pattern formation of microorganisms and blow-up phenomena in finite time. 
Key contributions on chemotactic collapse have been brought by authors in ora uni mum. Similarly, 
angiogenesis is the motion of endothelial cells to form new blood vessels from preexisting vessels to locally 
supply oxygen. During tumor growth, tumor cells secrete a set of substances to attract endothelial cells m 
□g. In angiogenesis models, chemotactic collapse was not observed nn ns]. In this study, the model equations 
(l|4| are used to describe chemotactic response of bacterial suspensions. The unstable agglomeration of 
bacterial cells at the surface leads to a descent of cell-rich plume patterns together with the formation of 
high-speed jets between counter-rotating vortices. 

Formation and stability of plumes result from the balance between chemotaxis, diffusion, and convection 
of bacteria. The particular impact of each mechanism still needs to be understood. Chemotaxis is known to 
bring instability in system and leads to aggregation, but may also have a stabilizing effect. 

The linear stability analysis of the chemotaxis-diffusion-convection system showed that a condition for 
linear instability depends on the taxis Rayleigh number Ra T [2]. The taxis Rayleigh number is defined 
as the ratio of buoyancy and viscosity forces times the ratio of momentum and cell diffusivity (table |2j) . 
Below a critical value ( Ra T < Ra c ), then the process remains stable. From experiments, several stages 
were observed starting from the upward bacteria accumulation and leading to hydrodynamic formation of 
plumes. A weakly nonlinear stability analysis was conducted to investigate the stability of hexagon and roll 
patterns formed by the system of equations (l|4) [17] . The hydrodynamic vortices formed by convection 
strengthen circulation of fluid and enhance the intake of oxygen into the solvent [T5] . Global existence for 
the chemotaxis-Stokes system for small initial bacterial population density was proved in m ■ Then, global 
existence for the chemotaxis-Navier-Stokes system for a large initial bacterial population density as well as 
global existence of 3D weak solutions for the chemotaxis-Stokes equations were proved in -2D3 ■ 

A detailed numerical study of the plume formation and merging that was related to the convergence of 
Rayleigh-Benard-type patterns was carried out in [2Tj . The shape and number of plumes can be controlled by 
initial bacterial population density. However, the sites of plumes have not been predicted. The convergence 
toward numerically stable stationary plumes for low and high initial density of cells is revealed in |2I| . 

In this paper, a computational model based on the finite element method is proposed aiming at investi- 
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Table 2: Representative dimensionless numbers involved in the double diffusive convection (DDC), 
chemotaxis-diffusion-convection (CDC), and Rayleigh-Benard convection (RBC). Subscripts - m , -t, -s, 

• o, ■b stand for mass, thermal, taxis, solute, oxygen, and bacterium, respectively; g is the acceleration due 
to gravity; /3t and /3 S are the coefficients for thermal and solute expansion, respectively; v is the kinematic 
viscosity. 



DDC 

CDC 

RBC 

Rayleigh 

number 

Ra T — 

Ut v 

D „ _5/3 s AsL 3 

D„ v 

D gV b nR(p b -p) L 3 

Ra r = 

D b p 

p g Pt AT L 3 

Ro, t = 

Dtp 

Prandtl 

number 


■i 

s 

II 

£1* 

FrT= ^- T 

Lewis 

number 

r Dr 

Ut - d7 

T D ° 

T D b 


Chemotaxis 

sensitivity 


q ^ dim ^air 

d 6 


Chemotaxis 

head 


|_| _ K ^0 L 2 

^air Db 



gating the behavior of the two-dimensional chemotaxis-diffusion-convection system. We present numerical 
examples of different states of the system: (i) diffusion dominant, (ii) chemotaxis dominant, and (iii) con¬ 
vection dominant with the formation of descending plumes. Furthermore, we show that the chemotaxis 
sensitivity (S) and the taxis Rayleigh number ( Ra T ) are two relevant parameters for the generation of in¬ 
stabilities. The effects of initial conditions and the initial cell density were also explored. Distinct initial 
settings lead to different solutions with plume patterns. Numerical tests show how the specified deterministic 
initial condition influences the overall behavior such as the number of plumes. In addition, the chemotaxis- 
diffusion-convection system is compared with the Rayleigh-Benard and double-diffuse convection systems. 
The dimensionless parameters introduced in previous papers were renamed for better readability (table [2]). 

In section [2] of the present paper, the mathematical formulation of the chemotaxis-diffusion-convection 
system is presented. Section [3] gives the numerical method. Section [4] demonstrates the role of chemotaxis, in 
addition to a discussion on the simulation results and comparison with other models of buoyant convection. 
Finally, section [5] presents some concluding remarks. 
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2 Mathematical model 

A mathematical model for the chemotaxis-diffusion-convection is proposed in [3] and reads as follows: 


du 

~dt 


+ u-Vuj + Vp - pV 2 u = -n V b g(pb - p) j, 

(i) 

V-u = 0, 

(2) 

Otl 

— + V- [un - D b Vn + S dim r(c)nVc] = 0, 

(3) 

dc 

— + V- (u c — DoVc ) = —n «r(c), 
ot 

(4) 


where u = (u,v) denotes the velocity field of water (solvent), p the pressure, p and p the water density and 
viscosity, n the areal number density of bacteria (i.e., number of bacteria per unit area in a 2D space), Vj, 
and pb the volume and volumetric mass density of a bacterium, c the concentration of oxygen, Vbg(pb — p )j 
the buoyancy force exerted by a bacterium on the fluid in the vertical direction (unit vector j), S dim the 
dimensional chemotaxis sensitivity, Db and Do the bacterium and oxygen diffusivities, n the bacterial oxygen 
consumption rate, and r(c) the dimensionless cut-off function for oxygen concentration. The cut-off function 
r(c) is defined by the step function 

1 if c > c*, /r s 

(5) 


r{c) = 


0 if c < c* 


where c* = 0.3. 

Bacteria are slightly denser than water and are diluted in the solvent, so that we consider (pb — p)/p <C 1 
and nVb <C 1, respectively. The consumption of oxygen is proportional to the bacterial population density 
n. Both n and c are advected by the fluid. When the oxygen concentration is lower than a threshold, the 
bacteria become quiescent [3], that is, they neither consume oxygen nor swim toward sites of higher oxygen 
concentrations. The dynamics of space filling, intercellular signaling, and quorum sensing are also neglected. 

The system of equations (jT|4]) with the boundary conditions introduced in previous papers (e.g. ®mm) 
is solved in a two-dimensional rectangular container (fl). The top boundary (T top ) represents the interface 
between liquid and air. On the free surface the concentration of oxygen is equal to the air concentration of 
oxygen (c a j r ) and the free tangential stress condition as well as the absence of bacterial flux are prescribed 
(figure [l]). Therefore: 


du 

— = 0, v = 0, c = c air , 
dy 


( 6 ) 


Sdimnr(c) Vc-n — Db Vn-n = 0 on T top , 

where n is the unit outward normal vector. On the container walls (fo), a no-slip boundary condition is 
prescribed and the fluxes of bacteria and oxygen equal zero: 


u = 0, v = 0, Vn-n = 0, Vc-n = 0 on T„ 


(7) 


A no-slip condition at the air-water interface would enable the formation of hydrodynamic instabilities. The 
effect of a moving boundary due to the advection caused by an incompressible fluid flow is to be explored. 
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Air 

C C a i r 

No flux for n 


Free surface condition: 

y-direction velocity set to 0 



Wall 


No flux for n, c 


No slip condition 


Figure 1: Boundary conditions for the system of equations ( l|4). The air-water interface, where the oxygen 
concentration is equal to that of air, is not crossed by bacteria; the fluid vertical velocity component equals 
zero and the fluid is assumed to be free of tangential stress. The container walls are impermeable to bacteria 
and oxygen; a no-slip condition is imposed. 


3 Computational model 


3.1 Scaling and setting for numerical simulations 

The characteristic length is defined by the container height h and the characteristic bacterial density by the 
average of the initial bacterial population density defined as 

no(x,t)dx. (8) 



This particular choice of characteristic bacterial density allows us to easily measure the total number of 
bacteria in each simulation for different initial distributions of bacteria. 

Dimensionless variables are defined as 


, = x , t 
X h' h 2 /D b ’ 
P 


n 




c = 


P = 


pD b /h 2 ’ D b /h' 


(9) 


Five dimensionless parameters given below characterize the hydrodynamic and chemotaxis transport equa¬ 
tions: 


V n g V b n 0 (p b -p)L :i SdimCa.h 

>=-=-, na T = ---. b=-—- 

D b Db /j, D b 

nnd L 2 D 0 

= -—, Le r =—— 

Cair D b D b 


( 10 ) 


where Pr T is the taxis Prandtl number, Ra T the taxis Rayleigh number (buoyancy-driven flow), Le T the taxis 
Lewis number, S the dimensionless chemotaxis sensitivity, and H the chemotaxis head. The taxis Prandtl, 
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Rayleigh, and Lewis numbers are analogous to the respective heat and mass Prandtl, Rayleigh, and Lewis 
numbers in heat and mass transfer (table [2]). The chemotaxis sensitivity (S) and head (H) characterize 
the chemotaxis system. Only Ra T and H depend on the characteristic length L and characteristic bacterial 
density nb- 

After removing the prime in dimensionless quantities, the set of dimensionless equations becomes 


dn 

dt 


d\i 

——h u-Vu — Pr T V 2 u + Pr T \7p = —Ra T Pr T n j, 
dt 

Vu = 0, 

+ u-Vn — V 2 n + S V-(r(c)nVc) = 0, 


( 11 ) 


dc 

dt 


+ u-Vc — Le T V 2 c = H n r(c). 


The system (18) is solved in a rectangular domain f l = [—£, i] x [0,1] with the initial conditions: 

u(x, 0) = u 0 ( x ) , n(x, 0) = n 0 (x), c(x, 0) = c 0 (x). 


( 12 ) 


On the top of the domain f1, the dimensionless boundary conditions are prescribed as 


Du 

— =0, f = 0, Sr (c) nVc-n — Vn-n = 0, c = 1, 
dy 

while on the other boundaries, we impose the dimensionless boundary conditions as follows 

u = 0, v = 0, Vn-n = 0, Vc-n = 0. 


(13) 


(14) 


In the following sections, we will refer to the hydrodynamic system for the first two equations in (18) and to 
the chemotaxis system for the last two equations in (18). 


3.2 Numerical methods 


The governing equations in (18) are solved using the finite element method. We adopt the biquadratic 


quadrilateral elements for the primitive variables u and the bilinear quadrilateral elements for the primitive 
variables p so as to satisfy the LBB (Ladyzhenskaya [22] - Babuska (23] - Brezzi [24]) stability condition. 
There are nine nodes in one biquadratic quadrilateral element and four nodes in one bilinear quadrilateral 
element. In each element, the function </> can be written as (f> = where 4>i are the nodal unknowns. 

To avoid the convective instability while solving the convection dominated flow equations, we adopt 
the idea of a streamline upwind/Petrov-Galerkin (SUPG) method 25- We implemented an inconsistent 
Petrov-Galerkin weighted residual scheme developed in [IS]. The resulting weak formulation reads as 


f dcj) dc/) d(j) dN d<j> dNd(j) 

P N i + + "TTyl + ' ^ & £ + T^iTy 


= Jm + H N f x 


0(j) 

+ 

dy 


V2 


)• 


(15) 


In the context of inconsistent formulation, we have two kinds of test functions. For the non-convective 
terms we choose the test function to be the shape function N. The test function W is only applied to the 
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convective term so as to add a numerical stabilizing term. The main purpose of employing a SUPG scheme 
is to add an amount of streamline artificial viscosity. The test function W is therefore rewritten in two parts 
W = N + B. B is called the biased part. On each node i, the biased part is Pj = TUj where j € {1,2} 
corresponds to the Cartesian coordinates and (iq, u 2 ) = (u,v). 

From the exact solution of the convection-diffusion equation in one dimension, following the derivation 
in 12 E], the constant r is determined as 


r = 


5 ( 7 i) u hi + S( 72) v h 2 
2 (u 2 + v 2 ) 


(16) 


We define 7 j = Uj 2 ] with (h\,h 2 ) = (Ax, Ay) being denoted as the grid sizes. Finally, the derived expression 
for 6( 7 ) is given below 


S(j) = 


2-cosh(7)-^ tanh(J)+i sinh(7) 

4 tanh(^)—sinh(7)—— sinh(7) tanh(^)’ 


cothm - i 


at end-nodes, 
at center-nodes. 


(17) 


For comparison purpose in section |4.4| the double diffusion system (26) and the Rayleigh Benard system 


(27) are solved using the software Freefem-)—(- S2z:. 
weak formulation of the problem. Taylor-Flood P 2 -I 
are used together with a characteristic/Galerkin formulation to stabilize the convection terms. 


The code uses a finite element method based on the 
1 elements are chosen in FreeFem-|—(-. These elements 


3.3 Numerical validation of the coupled NS-KS equations 

In this validation study, the following dimensionless differential equations accounting for the coupled 
Keller-Segel and incompressible viscous hydrodynamic equations are solved 


du 

at 


u-Vu — Pr r V 2 u + PryVp = —Ra T Pr T n j + f u , 
V-u = 0, 


dn 

at 


+ U'Vn — V 2 n + S V-(r(c)nVc) = /„ 


(18) 


dc 

at 


+ u-Vc — Le r V 2 c = H n r(c) + f c in fi. 


The physical properties are set at the constant values of Pr T = Ra T = S = r(c) = Le T = H = 1 in Cl = 
[0,1] x [0,1]. Equations in ( |18| are solved at At = 0.0001 in the continuously refined four meshes with Ace = 
Ay = 0.125,0.1,0.0625,0.05. The predicted errors between the simulated and exact solutions, which are 
u exact = — cos(7rx) sin(7T2/)e _27r *, v exact = sin(7rx) cos(7r?/)e _27r *, p = C\ — 0.25(cos(27rx) + cos(27ry))e _47r *, 
nexact = cos(7rx) cos(7 Ty)e~ 2 ™ 1 and c exac t = cos(7rx) cos(7t y)e~ 27T *, are cast in their L 2 — norms. The source 
terms f u , / ra , and f c are derived from the previous exact solutions. From the predicted error norms, the 
spatial rates of convergence are plotted in Fig. [2] 

The good agreement between the exact and simulated results and states of convergence demonstrate the 
applicability of the proposed SUPG scheme and the flow solver described in section |3.2| to investigate the 
chemotactic phenomena in hydrodynamic environment. 









(a) 



(b) 


Figure 2: The computed rates of convergence (roc) for the coupled set of Navier-Stokes and Keller-Segel 
equations, (a) roc= 3.95382 for u, roc= 3.92594 for v, roc= 1.96484 for p; (b) roc= 2.81384 for n, roc= 
1.80646 for c. 9 















4 Numerical results and discussion 


The linear stability analysis of the system (18) showed that the steady state becomes unstable for a range of 
physical parameters [2]. For sufficiently large characteristic bacterial density tTq and Rayleigh number Ra T , 
hydrodynamic instabilities appear in the region near the surface at which the bacterial density is high. This 
instability may be related to the Rayleigh-Benard instability occurring in thermal convection [3|, [T7] - This 
instability develops into a descending family of bacterium-rich plumes and leads possibly to the formation 
of convection cells (figure [3] (a)). 

When no is small, small perturbations of the velocity field are damped due to stabilizing effects of 
viscous friction and chemotaxis and therefore convective motion is negligible. It is shown that for low initial 
bacterium cell density nd the system (18) evolves to a steady and homogeneous state in the horizontal 
direction governed by the chemotaxis system EQI2j. Like in the Keller-Segel or angiogenesis system, the 
random diffusion of cells in this case is balanced by the chemotaxis sensitivity of cells. 


4.1 Descending plumes 

In the present numerical simulations, the descending plumes can be described using three phases. In the 
first phase (figure [4] (a)-(c)), chemotaxis is a dominant mechanism. As bacteria consume oxygen, an oxygen 
concentration gradient is created that in turn provokes a bacterium motion toward the open surface where 
the oxygen is abundant. Bacteria chemotaxis causes the fluid to set in motion and counter rotating vortices 
to form (figure[3](a)). Quickly, the bacterial density n becomes quasi-homogeneous in the horizontal direction 
and is structured in layers in the vertical direction. A three-layer configuration is induced. A layer of higher 
concentration of bacteria forms below the air surface: the stack layer. When bacteria have migrated in the 
stack layer, a depletion layer is generated. At the container bottom, some bacteria are inactive, because the 
oxygen concentration decreases below a certain threshold. An inactive layer is established. 

The second phase (figure [4] (d)-(e)) exhibits high bacterial density on the surface. As bacteria swim 
toward air-supply region, the bacterial density in the stack layer increases. Consequently, advection from 
the counter rotating vortices becomes significant and brings in perturbations to the stack layer (figure [5| . 
Therefore, advection causes instabilities to occur in the bacterial density distribution in the stack layer that 
sets the fluid in motion. At the same time, the oxygen concentration falls to its minimum at the bottom. 

In the third phase (figure [4] (f)-(h)), in fluid regions with a greater bacterial density in the stack layer, 
buoyancy force constrains bacteria to descend in the fluid. As a result, descending plumes of bacteria develop 
at these particular locations. This mechanism is analogous to the Rayleigh-Benard instability in heat transfer 
problems where the fluid with a higher temperature and thus a lighter density than that above it rises to 
develop plumes of hot fluid. 


4.2 Stabilizing effect of chemotaxis 


When the average initial cell density no is large, hydrodynamic instabilities appear in the system (18). This 
section is aimed at estimating nd and dimensionless governing parameters for hydrodynamic instabilities 
to appear as well as at analyzing the time scales for each of the three competitive physical mechanisms: 
(1) chemotaxis, (2) diffusion, and (3) convection of bacteria. 

Convection is determined by the properties of the hydrodynamic system as well as the container height 
and nd via the taxis Rayleigh number Ra T . When the Rayleigh number is below the critical value for the 
hydrodynamic system of interest, motion of bacteria is primarily governed by diffusion and chemotaxis; when 
the Rayleigh number exceeds a critical value, bacterial taxis is primarily governed by convection [2]. 
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Figure 4: Evolution of the cell density n at different times. Descending plumes of bacteria develop from an 
initial randomly distributed bacterium population. 


Buoyancy force enables a fluid volume with a low bacterial density to ascend and a fluid volume with a 
high bacterial density to descend. Nevertheless, bacterium chemotaxis and friction dampen altogether the 
displacement. 

The time scale for bacterium diffusion over the length scale h is given by 

T diff := (19) 

The buoyancy force is balanced by friction in the fluid. Therefore, the time scale for the convective displace¬ 
ment of bacteria over the length scale h can be defined by 


T 


conv 


_a*_ 

ghnoV b (p b - p)' 


( 20 ) 


The chemotaxis system is controlled by a competition between diffusion and chemotaxis of bacteria. There¬ 
fore, the chemotaxis time scale is expressed as 


= Db 
S dim bv TIq 


( 21 ) 


In a 
time 


convection-dominant process, the convection time scale is smaller than the diffusion and chemotaxis 
scales (table [3]). Hence: 


g h 3 n 0 Vb(pb-p) _ D ^ , 

-—- = tia T > I. 

D b p 


( 22 ) 
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Figure 5: Evolution of the cell density number n at the surface. In the initial stage of the chemotaxis- 
diffusion-convection, n is homogeneous in the horizontal direction. As the cell density n increases, hydrody¬ 
namic instability arises and the bacterial density becomes higher at some particular points. 
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Table 3: Phenomenological analysis based on time scales of the three competitive mechanisms: chemotaxis, 
diffusion, and convection of bacteria. Both oxygen and bacterial diffusion favor fluid homogenization. 


Dominant 

convection 

Dominant 

diffusion 

Dominant 

aerotaxis 

T C onv ^ Tdiff 

Tdiff < T t 

T T < Tdiff 

and Tconv ^ T r 5 

SH < 1 

SH > 1 

i.e., Ra T > 1 



and Ra T /SH > 1 




In particular, Ra T > Ra c . The value of the critical Rayleigh number associated with the convection Ra c 
given by j2] in their linear stability analysis is described by the solution of an ordinary differential system. 
The second condition on the time scales leads to the following inequality 


ghD b V b (p b - p) _ Ra T 
S dim ftp S H 


> 1 . 


(23) 


Similarly to the critical Rayleigh number ( Ra c ), a critical number (SH) C can be introduced. 

On the other hand, the chemotactic motion is predominant when the chemotaxis time scale is smaller 
than the diffusive time scale and convection is negligible: 


S dim K TIq h 

D b 2 


= SH > 1. 


(24) 


The effects of Ra c and the product S H were tested subject to the random initial condition. The results 
are shown in figure [6] for Pr T = 500 and Le T = 5. Our numerical results agree with the linear stability 
analysis carried out in [2]. Ra c at first falls as S H increases to reach its minimum value and then rises 
again. A sufficient decrease or increase of S H may promote stabilization as both taxis and diffusion operate 
as the fluid-homogeneization factors, either directly (S H < 1) or indirectly (S H > 1). When S H is small, 
stabilizing effect is due to bacterial diffusion and the solution is of the diffusive type. When S H is large, 
stabilizing effect results from the bacterial taxis and the solution is of the chemotactic type. 

The taxis Rayleigh number ( Ra T ) that characterizes the competition between diffusion and convection 
plays the same role as the Rayleigh number plays in classical convection. When the Rayleigh number 
increases, the gravitational force becomes predominant. When S H rises, competition between chemotaxis 
and convection of bacteria becomes stronger. The condition (231 suggests that Ra c increases linearly with 


respect to SH, as illustrated in figure [bj We can affirm that chemotaxis has a stabilizing effect on the 
differential system of current interest. 


4.3 Distribution and number of plumes and initial conditions 

4.3.1 Position and spacing of plumes 

The exact localization of plume generation was investigated. Under the random initial condition introduced 
in EB, the location is hard to be predicted. We consider deterministic initial conditions to investigate the 
formation and location of descending plumes. Many initial conditions among the set of tested distribution 
of bacterial settings are able to trigger convective patterns. 

All numerical results are computed at Le T = 5, Pr T = 500, Ra T = 2000, S = 10, and H = 4 to ensure the 
formation of bacterium rich plumes. We first consider a profile given by the form of wave function cos(27ra;/3) 
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SH 1 4 


Figure 6: The stable and unstable regions of the system (18) are plotted in terms of the critical taxis 
Rayleigh number Ra c and the product of the dimensionless chemotaxis sensitivity and chemotaxis head S H. 
The points correspond to the predicted values of Ra c . The red line corresponds to a fit of data in two parts: 
on the left side it is fitted by a power function and on the right side by a linear function. In the unstable 
region, the predicted solution n is of the convective type (b). In the stable region the predicted solution n 
can be of the diffusive type (a) or the chemotactic type (c). 
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(C) (d) 


Figure 7: Numerical results for n with respect to the deterministic initial conditions. At the initial time, bac¬ 
terial density is higher in the upper layer, (a) initial condition is no = [ g’ 5 a ^<o^+oTcos( 27 rx 3 / 3 )> 0 3 ) simulation 

result at t = 1.2 for the initial condition given in (a) ; (c) initial condition is no = [ g’ 5 5-0 i°cos( 2 rra/ 3 ) > 

(d) simulation result at t = 1.2 for the initial condition given in (c). 


to determine the two initial layers of bacteria. Simulation results and initial conditions are shown in figure 
[T] The upper layer has a higher bacterial density than the lower layer. The plumes form at the initial 
location of the crest of the wave function where there is a larger number of bacteria. By reversing the initial 
conditions with a denser inferior layer (figure [8]), the plumes also form at the same location. We can observe 
that distinct initial conditions can lead to the formation of a very similar pattern (figures [7]|sj). 

We vary the profile of the curve between the upper and lower layers in the initial condition with different 
wave functions cos(37ra;/2), cos(57ra/2), and cos(27ra;). The wavenumber is defined as the spatial frequency 
of waves per unit distance. The wavelength is defined as the domain horizontal length divided by the 
wavenumber. Increasing the wavenumber of the initial wave function causes the formation of additional 
plumes to occur (figures [9] and 10). In all simualtions, the locations of the plumes correspond to sites of 
initially greater local bacterial density. 

When the wavenumber in the initial profile is large enough, the number of plumes is not equal to the 
wavenumber fixed by the initial condition. Three plumes, at most, form (figures 11 to 14). Nonetheless, 


more than three plumes can form, but they merge later (figure [l3| ). Plume merging was previously observed 
in [51]. The plume merging mechanism is analogous to that of the Rayleigh-Benard flow. Spacing between 
plumes seems to be intrinsic to the system and the position can only be predicted in very specific cases such 
as the examples given above. 

The randomly perturbed bacterial density is defined as follows |21| : 


n(x) = 0.8 + 0.2e(x), 


(25) 


with e being a random number with a uniform probability distribution over [0,1]. Simulation results with 
the randomly perturbed initial condition are given in figure [15] We note that the plume-to-plume spacing 
is not fixed but the number of plumes in the solution remains the same as that with the solution obtained 
subject to the deterministic initial condition illustrated in figures 

In figures |7| |15[ we can see plumes forming at the border of the domain. These border plumes seem to be 
caused by the no-slip boundary condition on the velocity depending on the arrangement of the convection 
cells. Bacteria first agglomerate at the surface near the wall. Because the fluid velocity is equal to zero at 
the wall boundary, a large amount of bacteria remains close to the wall. At a location away from the wall, 


16 
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Figure 8: Numerical results for n with respect to the deterministic initial conditions. At the initial time, bac¬ 
terial density is higher in the upper layer, (a) initial condition is no = [ 1 5 ’j / <o' 5 +o _ i cos° 27 ^c/ 3 ) 3 ^ ’ (b) simulation 

result at t = 1.2 for the initial condition given in (a) ; (c) initial condition is n 0 = [ 5-0 icos( 2 ir®/ 3 ) 3 ^ 

(d) simulation result at t = 1.2 for the initial condition given in (c). 



(c) (d) 


Figure 9: Numerical results for n with respect to the deterministic initial conditions. At the initial time, bac¬ 
terial density is higher in the upper layer, (a) initial condition is no = [ g’ 5 ^<o~ 5 +o i°cos( 57 ra/ 2 ) > 0 3 ) simulation 

result at t = 1.2 for the initial condition given in (a) ; (c) initial condition is no = [ g ’ 5 V ~y K g 5^0 i°cos(57ra/2) > 
(d) simulation result at t = 1.2 for the initial condition given in (c). 
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Figure 10: Numerical results for n with respect to the deterministic initial conditions. At the initial time, bac¬ 
terial density is higher in the lower layer, (a) initial condition is no = [ 0 3 ) simulation 

result at t = 1.2 for the initial condition given in (a) ; (c) initial condition is n 0 = [ 5-0 icos°5w®/2) 2 ^> 

(d) simulation result at t = 1.2 for the initial condition given in (c). 



(c) (d) 


Figure 11: Numerical results for n with respect to the deterministic initial conditions. At the initial time, bac¬ 
terial density is higher in the upper layer, (a) initial condition is no = [ q’ 5 ^< 0^+0 i°cos( 37 rx 2 / 2 ) > 0 3 ) simulation 

result at t = 1.2 for the initial condition given in (a) ; (c) initial condition is no = [ q ’ 5 5^0 i°cos( 37 ra/ 2 ) ’ 

(d) simulation result at t = 1.2 for the initial condition given in (c). 
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Figure 12: Numerical results for n with respect to the deterministic initial conditions. At the initial time, bac¬ 
terial density is higher in the lower layer, (a) initial condition is no = [ ' 5 ’j / <^ 5 +o'i'cos( 37 fJ/ 2 ) 2 ^ 0 3 ) simulation 

result at t = 1.2 for the initial condition given in (a) ; (c) initial condition is n 0 = [ 5-0 icos°3if®/2) 2 ^ 

(d) simulation result at t = 1.2 for the initial condition given in (c). 



(c) (d) 


Figure 13: Numerical results for n with respect to the deterministic initial conditions. At the initial time, 
bacterial density is higher in the upper layer, (a) initial condition is tiq = [ q’ 5 !/ ^<o~ 5 +o i°cos( 2 rra;) ’ (b) simula¬ 
tion result at t = 1.2 for the initial condition given in (a) ; (c) initial condition is no = [ q ' 5 ^<0 5^0 i°cos( 2 ff®)> 
(d) simulation result at t = 1.2 for the initial condition given in (c). 
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Figure 14: Numerical results for n with respect to the deterministic initial conditions. At the initial time, 
bacterial density is higher in the lower layer, (a) initial condition is no = [ i 5 ’y<^ 5 +o"icos° 27 ra) X ^ (b) simula¬ 
tion result at t = 1.2 for the initial condition given in (a) ; (c) initial condition is no — [ i cos^ 27 ra) X ^> 

(d) simulation result at t = 1.2 for the initial condition given in (c). 



Figure 15: Simulations subject to the initial condition given in (25) (left) and the corresponding numerical 
results for n at t = 1.2 (right). 
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Table 4: The predicted number of plumes, wavenumber and wavelength for £ = 2. Each row corresponds to 

lly perturbed initial conditk 
wavelength wavenumber 


a different simulation result subject to the randomly perturbed initial condition given in (25). 

number of plumes 


3.14 

3.14 

3.14 

3.14 


Table 5: The predicted number of plumes, wavenumber and wavelength for £ = 3. Each simulation result is 
subject to the random initial condition given in (251. 

number of plumes wavelength wavenumber 


3 2 3.14 

3 2 3.14 

4 1.5 4.19 

3 2 3.14 


the velocity is non-zero, thus bacteria descend into the fluid and form a plume at the border. 

The location of the plumes may only be predicted in very simple tests for which the wavelengths of 
initial conditions are lower than the wavelength of the process. Despite the difficulty to predict the site of 
plume formation, the system presents a dominant wavelength of the fingering instability. The wavelength 
and wavenumber are redefined such as those given in |28j . The wavelength is the length of the domain 
divided by the number of plumes and the wavenumber 2n divided by the wavelength. The wavelength and 
wavenumber are similar for all simulation tests whatever the length of the domain is (tables [4] to [7| . 

4.3.2 Growth of plumes 

We now arbitrarily define the plumes by the isoline n = 1. This choice allows us to track the stack layer 
and the plumes being formed. The layer where n > 1 represents the layer with a high bacterial density from 
which plumes form and descend in the fluid. On the other hand, the layer with n < 1 corresponds to the 
depletion layer, from which bacteria have migrated to the stack layer. 

We define the growth rate of the plume amplitudes by G = ( A t — A t _i)/At. The plume amplitude A t 
is computed by measuring the distance from the surface to the tip of the plume at time t (table [8]) and At 
is the time increment. Data obtained can be interpolated by the exponential function g(t) = expjcrf + /3} 
as the amplitudes of the descending plumes grow exponentially (figure |16[ ) . The exponential growth is well 
understood in the conventional Rayleigh-Taylor convection. Subsequent to the exponential growth phase, 


Table 6: The predicted number of plumes, wavenumber and wavelength for £ = 4. Each simulation result is 
subject to the random initial condition given in (251. 

number of plumes wavelength wavenumber 


5 

1.6 

3.93 

5 

1.6 

3.93 

5 

1.6 

3.93 

5 

1.6 

3.93 
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Table 7: The predicted number of plumes, wavenumber and wavelength for £ = 5. Each simulation result is 
subject to the random initial condition given in (251. 

number of plumes wavelength wavenumber 


2 

1.67 

2 

2 


3.14 

3.77 

3.14 

3.14 


Table 8: The predicted growth rate G of the plume amplitudes A t for two simulations subject to the random 

(a) (b) 


initial condition (251. 


At 

t 

G 

At 

t 

G 

0.18 

0.16 


0.18 

0.16 


0.19 

0.18 

0.46 

0.20 

0.18 

1.07 

0.20 

0.19 

0.76 

0.21 

0.19 

1.07 

0.21 

0.20 

0.92 

0.23 

0.20 

1.99 

0.23 

0.21 

1.53 

0.28 

0.21 

3.53 

0.26 

0.22 

2.91 

0.37 

0.22 

8.13 

0.34 

0.24 

6.44 

0.57 

0.24 

16.34 

0.51 

0.25 

13.81 

0.78 

0.25 

17.12 

0.72 

0.26 

17.81 

0.87 

0.26 

8.13 

0.83 

0.27 

9.67 

0.93 

0.27 

4.91 


the growth rate of amplitude decreases, as plumes get closer to the bottom of the container. 


4.4 Comparison with other buoyancy-driven convections 

The chemotaxis-diffusion-convection system has many features similar to other well known buoyancy-driven 
flows, such as the double diffusive and Rayleigh Benard convection. 

Double diffusive convection occurs in a fluid containing at least two components with different diffusivities. 
A destabilizing component diffuses faster than the stabilizing one [S!J|. The distinct diffusivities yield a density 
difference capable of driving the motion of fluid [30l . 

Comparison with a classical example of double diffusive phenomenon in oceanography can be carried out 
by considering two superposed fluid layers with a specific combination of temperature (T) and the solute 
concentration (s), namely the salinity. The system of dimensionless equations of the double diffuse problem 
is the following: 

(9u 

—+ (u-V)u-Pr T V 2 u+Pr T Vp = -Pr T (Ra m s-Ra T T)j , 
at 

V-u = 0, 

dT 2 (26) 

—+ u-vp- v 2 r = o, 

at 

ds 

— + u-Vs — Ler V 2 s = 0. 
at 
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Figure 16: Plot of the growth rate G of plume amplitude (dots) corresponding to the data in 
growth rate is interpolated by the exponential function g(t) = e at+l3 (line). 
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table [8] The 








where Rax and Ra m are the thermal and mass Rayleigh number, respectively, Prx the Prandtl number, and 
Lex the Lewis number (table [2]). 

The stability of the system that exhibits a diffusive and a finger regime depends on both Rayleigh number 
types. In the finger regime, a small perturbation at the interface between the layers develops into a pattern 
of descending fingers, for instance salt fingers. 

Another buoyancy-driven convection is the Rayleigh-Benard convection that arises by fluid in a reservoir 
heated from below. Convection results from thermal gradient. The set of equations is as follows: 


—+(u-V)u— Prx \7 2 u+Prx Vp = —Prx RaxT j, 

Vu = 0, 

QT 

— + U-VT - V 2 T = 0, 
at 

T = 1 at bottom, T = 0 at top. 


(27) 


Rax and Prx are the Rayleigh and Prandtl number, respectively. The balance between the gravitational 
and viscous forces is expressed by the Rayleigh number Rax ■ When Rax is larger than a critical value that 
can be obtained analytically, convective patterns appear m- 

The double diffusive convection equations in (26) and the Rayleigh-Benard convection equations in(27) 
are solved using our numerical method (figures 17 and 18). Plumes formed in the Rayleigh-Benard convection 
process are similar in shape to plumes of the chemotaxis-diffusion-convection. However, double diffusive and 
Rayleigh-Benard convections exhibit both ascending and descending plumes, whereas chemotaxis-diffusion- 
convection is only characterized by descending plumes. 

In each problem, the effective Rayleigh numbers depend on the temperature difference between the 
opposite fluid domain surfaces, the gradients of temperature and salinity between the fluid layers, and 
the difference of density between the bacteria and solvent, for the Rayleigh-Benard, double diffusive, and 
chemotaxis-diffusion-convection systems, respectively. All of the three diffusion-convection processes have 
similar and distinct features (table [9]) . The involved dimensionless parameters are listed in table [2j 

Figure [3] shows that the arrangement of the convection cell structure can be the same for the three systems 
mentioned above. A particular position of the domain will be in a clockwise cell in some simulations and 
counter-clockwise cell in others. This behavior is also observed in Rayleigh-Benard convection simulations. 

In the chemotaxis-diffusion-convection system, chemotaxis plays an essential role in the early stage, as 
it organizes the fluid domain. In a rectangular domain, starting from a given initial condition, aerotaxis of 
bacteria builds quasi-homogeneous layers in the horizontal direction, similarly to the double diffusive m- 
The multi-layered fluid creates a density gradient between the top of the stack layer and the bottom of the 
depletion layer that is similar to the temperature gradient set by the Dirichlet boundary condition imposed in 
the Rayleigh-Benard case. The chemotaxis-diffusion-convection system evolves itself to a proper condition 
that leads to instabilities. 

Chemotaxis, that is not present in the other systems, brings flexibility in the choice of initial as well as 
boundary conditions. An analogy between the differential systems should stimulate further mathematical 
analysis of the system and a better understanding of the role of chemotaxis in convection problems. 


5 Concluding remarks 

In this paper, we have studied the chemotaxis-diffusion-convection system with the focus on the differential 
system rather than on the experimental settings. Our studies are based on the numerical implementation 
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(a) 



(b) 


Figure 17: Plot of the solution s in the double diffusive system (261. Red area designates the region rich in 
s, while yellow area is poor in s. Descending plumes rich in s and ascending plumes poor in s form at the 
interface of two layers of fluid. The green color designates a region where the concentration of s is slightly 
higher than that in the yellow area due to diffusion, (a) Initial condition: s(x,y,t = 0) = 2/7, if y < 0.5, 
and s(x, y) = 1.0, if y > 0.5. (b) s at time t= 0.1 for Rqt = 8000, Ra m = Ra T /6. 2, Prr = 7 and Ler = 0.01. 



Figure 18: Plot of the solution T in the Rayleigh-Benard system (271. Descending plumes of lower temper¬ 
ature and ascending plumes of higher temperature form in the fluid. Rax = 6000, Prx = 50. 
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Table 9: Recapitulative of the physical mechanisms involved in the double diffusive convection (DDC), 
chemotaxis-diffusion-convection (CDC) and the Rayleigh-Benard convection (RBC). 



DDC 

CDC 

RBC 

Hydrodynamics 

buoyancy (pVg) 

V-u = 0 

7j^-=PrxV 2 u— PrrVp 
—Pr T ( Ra m s-Ra T T )j 

V-u = 0 

j^-=Pr r V 2 u — Pr T Vp 
—Ra T Pr T n j 

V-u = 0 

-^-=PrT V 2 u —Ptt Vp 
—Pvt Rut j 

Diffusion 

DT 

—— = V 2 T 

D t 

Ds 

—- = LerV 2 s 

D t 

Dn 

— = V 2 n 

Df 

DT 

— = V 2 T 

D t 

Chemotaxis 

0 

1~)77 

— = -SV-(nVc) 

Df v ' 

Dc ,, 

— = Le r V c — Hn 

Df 

0 

Convection 

DT 

— = u-vr 

D t 

Ds 

— = u-Vs 

Df 

Dn 

= u-Vn 

Df 

Dn 

■=— = u-Vn 

Df 

I.C. for fingers 

Layers 

Any 

Layers 

B.C. for mass 

Neumann 

Neumann 

Dirichlet 

Physical regime 

Diffusive & 

convective 

Diffusive, 

convective & 

chemotactic 

Diffusive & 

convective 
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of the upwind finite element method with inconsistent Petrov-Galerkin weighted residual scheme to solve 
the coupled convective chemotaxis-fluid equations. Several simulation examples have been presented. They 
exhibit physically different spatially organized convection patterns. 

The chemotaxis-diffusion-convection system can be described by three regimes. In the convective regime, 
the formation of plume patterns proceeds a temporal sequence of development stages. First, a spatial layered 
structure is created. Bacteria agglomerate in the upper stack layer, thus forming a depletion layer, where 
the bacterial density is very low. Subsequently, bacterial convection strengthens with time and instabilities 
in the stack layer appear. Finally, plumes of bacterial falling in the fluid are generated. The growth rate 
of the plumes is of the exponential type. The wavelength of the growing instabilities is defined from the 
parameters of the problem. Initial conditions appear to have a small influence on the plume number. The 
location of plumes can only be predicted when a very simple initial condition is set up. 

In the diffusive and chemotactic regimes, the chemotaxis system has a stabilizing effect on the fluid. When 
the taxis Rayleigh number increases, the gravitational force becomes dominant and instabilities occur and 
grow. From the phenomenological analysis and the numerical results, the critical Rayleigh number increases 
proportionally with the product of the chemotaxis sensitivity and head, when the chemotaxis sensitivity is 
large. 

In all buoyancy-driven flows (chemotaxis-diffusion, double diffusive, and Rayleigh-Benard convection), 
the dimensionless differential systems and the convection patterns are similar. Analogy between these types 
of convection should launch further analytical studies of the chemotaxis-diffusion-convection problem. 
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